*************************************************************
* The following code is part of the replication archive of  * 
* Genovese, Kern & Martin (2016) "Policy Alternation", ISQ  * 
*************************************************************

/* The code replicates the analyses that are based on an MSTAR model (the model originally used in the Ward and Cao CPS paper). 
 The estimates show the effect of spatial lag of energy/environmental subsidies on green taxes.  
 Note: The Stata MSTAR code below is discussed in Franzese, Hays and Kachi (2009). */

***************** 
* Summary stats *
*****************

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies.dta", clear
sutex greentax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp subsidy_rgdppc , nobs minmax   /* sl_gled_tax sl_gled_sub sl_dtrade_tax sl_dtrade_sub sl_igo_tax sl_igo_sub */

************
* Figure 2 *
************

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies.dta", clear
cd "/Users/.../gkm_replication_isq/green taxes and subsidies/Figures"
graph set window fontface "Times New Roman"

corr greentax_pc_cons subsidy if country=="Poland" 
corr greentax_pc_cons subsidy if country=="Czech"
corr greentax_pc_cons subsidy if country=="Austria" 
corr greentax_pc_cons subsidy if country=="Germany"

#delimit ;
twoway (line greentax_pc_cons year, yaxis(1)   bcol(gs12))  (line subsidy_rgdppc year,  yaxis(2)  lpattern(dash) lcol(red)) if country=="Poland" |country=="Czech"  ,  by(country, row(1) plotregion(fcolor(white)) graphregion(fcolor(white)) note(" ")  r2title("Level of Green Subsidy per capita"))   
ylab(0(300)960)
ytitle("Level of Green Tax per capita")
ylab(0(2000)9000,  axis(2))
xlab(1994 "94"  1996 "96"  1998 "98" 2000 "00" 2002 "02" 2004 "04")
legend(label(1 "Level of Green Tax") label(2 "Level of Green Subsidy"))
legend(cols(3) size(med) symysize(3) span ) legend(region(lwidth(none)))
plotregion(fcolor(white)) graphregion(fcolor(white)) saving(pol_czech, replace) ;
#delimit cr
 
#delimit ;
twoway (line greentax_pc_cons year, yaxis(1)   bcol(gs12))  (line subsidy_rgdppc year,  yaxis(2)  lpattern(dash) lcol(red)) if country=="Austria" |country=="Germany"   ,  by(country, row(1) plotregion(fcolor(white)) graphregion(fcolor(white)) note(" ") r2title("Level of Green Subsidy per capita"))   
ylab(0(300)960)
ytitle("Level of Green Tax per capita")
ylab(0(2000)9000,  axis(2))
xlab(1994 "94"  1996 "96"  1998 "98" 2000 "00" 2002 "02" 2004 "04")
legend(label(1 "Level of Green Tax") label(2 "Level of Green Subsidy"))
legend(cols(3) size(med) symysize(3) span ) legend(region(lwidth(none)))
plotregion(fcolor(white)) graphregion(fcolor(none)) saving(aus_ger, replace)  ;
#delimit cr

gr combine  aus_ger.gph pol_czech.gph, xcommon  plotregion(fcolor(white)) row(2)  graphregion(fcolor(white)) iscale(.55)  saving(Figure2, replace)

*****************************************
* Close countries:  Results in Table 1  *
*****************************************

* MSTAR model with 3 spatial matrices: geographical distance, IGOs and  trade dyads + subsidy 

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 sigma
tempvar A rSL1 rSL2 rSL3
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1 - p2W2 - p3W3
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3'-`mu', 0, `sigma'))
end

**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 

**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc   

mkmat $Y, matrix(Y)

matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X) (rho1:) (rho2:) (rho3:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo D1

drop SL11
drop SL21
drop SL31

********************************************

* MSTAR model with spatially lagged subsidy (close countries)

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 rho4 sigma
tempvar A rSL1 rSL2 rSL3 rSL4
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
gen double `rSL4'=`rho3'*SL4
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
scalar p4 = `rho4'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix p4W4 = p4*W4
matrix IpW = I_n - p1W1 - p2W2 - p3W3 -  p4W4
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3' - `rSL4'-`mu', 0, `sigma'))
end


**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag_close.dta", name(W4) standardize 

**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc  
global Z l_subsidy_rgdppc

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)

matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)
matrix SL4 = W4*Z
svmat SL4, n(SL4)

spatgsa  greentax_pc_cons , weights(W1) moran /* positive/negative coeff: nearby regions tend to exhibit similar/dissimilar values of Y*/
spatgsa  l_subsidy_rgdppc , weights(W4) moran /* positive/negative coeff: nearby regions tend to exhibit similar/dissimilar values of Y*/

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X) (rho1:) (rho2:) (rho3:) (rho4:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo D2

drop SL11
drop SL21
drop SL31

********************************************

* MSTAR model with interaction with spatially lagged subsidyXeconomic integration (flows)

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1 rho2 rho3 rho4 sigma
tempvar A rSL1 rSL2 rSL3 rSL4
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
gen double `rSL3'=`rho3'*SL3
gen double `rSL4'=`rho4'*SL4
scalar p1 = `rho1'
scalar p2 = `rho2'
scalar p3 = `rho3'
scalar p4 = `rho4'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix p3W3 = p3*W3
matrix p4W4 = p4*W4
matrix IpW = I_n - p1W1 - p2W2 - p3W3 -  p4W4
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'- `rSL3' - `rSL4'-`mu', 0, `sigma'))
end


**********************
*Open Data For Weights
**********************

spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag.dta", name(W1) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/trade_in2_lag.dta", name(W2) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/IGO_en_in2_lag.dta", name(W3) standardize 
spatwmat using "/Users/.../gkm_replication_isq/green taxes and subsidies/spatial weights/gled_dist_in2_lag_close.dta", name(W4) standardize 

**************************
*Open Data for Regression
**************************

drop _all

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear 

global nobs = 233 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y greentax_pc_cons
global X lag_gtax_pc_cons  rgdppc rgdppcsq unemp unempsq i_itax_per_cap actual_economic_flows rile_const_w env_const_w env_const_w_sq green_party en_prod_rgdp l_subsidy_rgdppc  
global Z l_subsidy_rgdppc

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)


matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W2*Y
svmat SL2, n(SL2)
matrix SL3 = W3*Y
svmat SL3, n(SL3)
matrix SL4 = W4*Z
svmat SL4, n(SL4)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X c.SL4#c.actual_economic_flows) (rho1:) (rho2:) (rho3:) (rho4:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo D3

************
* Figure 3 *
************

use "/Users/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies lag.dta", clear

reg greentax_pc_cons   lag_gtax_pc_cons   l_subsidy_rgdppc   rgdppc rgdppcsq  unemp  unempsq  i_itax_per_cap  actual_economic_flows  rile_const_w env_const_w   env_const_w_sq  green_party en_prod_rgdp     sl_gled_tax   sl_igo_tax  sl_dtrade_tax c.l_sl_neigh_subsidy_rgdppc##c.actual_economic_flows 
margins, at(l_sl_neigh_subsidy_rgdppc =(0(50)8000) actual_economic_flows =(35)) post
estimates store marginplot1

reg greentax_pc_cons   lag_gtax_pc_cons   l_subsidy_rgdppc   rgdppc rgdppcsq  unemp  unempsq  i_itax_per_cap  actual_economic_flows  rile_const_w env_const_w   env_const_w_sq  green_party en_prod_rgdp     sl_gled_tax   sl_igo_tax  sl_dtrade_tax c.l_sl_neigh_subsidy_rgdppc##c.actual_economic_flows 
margins, at(l_sl_neigh_subsidy_rgdppc =(0(50)8000) actual_economic_flows =(95)) post
estimates store marginplot2

cd "/Users/.../gkm_replication_isq/green taxes and subsidies/Figures"

graph set window fontface "Times New Roman"

coefplot  (marginplot1,  lpattern(ldash) lcolor(gs13))  , at ytitle(Marginal Effect of Economic Flows) xtitle(Spatially Weighted Green Subsidies) title("Low Levels of Economic Flows") recast(line) lwidth(*2) ciopts(recast(rline) lpattern(dash)) ylabel(500(100)800,nogrid) xsca(alt) levels(90) yline(0, lcolor(red)) plotregion(lcolor(black)) scheme(s2mono) graphregion(fcolor(white))   saving(marginplot1, replace)
coefplot  (marginplot2,  lpattern(ldash) lcolor(gs13))  , at ytitle(Marginal Effect of Economic Flows) xtitle(Spatially Weighted Green Subsidies) title("High Levels of Economic Flows") recast(line) lwidth(*2) ciopts(recast(rline) lpattern(dash)) ylabel(500(100)800,nogrid) xsca(alt) levels(90) yline(0, lcolor(red)) plotregion(lcolor(black)) scheme(s2mono) graphregion(fcolor(white)) saving(marginplot2, replace) 

twoway histogram  l_sl_neigh_subsidy_rgdppc if actual_economic_flows<76.5, freq gap(5) color(grey14) ylabel(0(25)50,nogrid) plotregion(lcolor(black)) xlabel(0(2000)8000,nogrid) fysize(22)  xtitle(Spatially Weighted Green Subsidies)  graphregion(fcolor(white)) saving(hist1, replace) 
twoway histogram  l_sl_neigh_subsidy_rgdppc if actual_economic_flows>76.5, freq gap(5) color(grey14)  ylabel(0(25)50,nogrid) plotregion(lcolor(black))  xlabel(0(2000)8000,nogrid) fysize(22) xtitle(Spatially Weighted Green Subsidies)  graphregion(fcolor(white))  saving(hist2, replace) 

graph combine marginplot1.gph hist1.gph, hole(2 4) imargin(0 0 1 0) graphregion(margin(l=10 r=10)) graphregion(fcolor(white))  saving(pic1, replace)
graph combine marginplot2.gph hist2.gph, hole(2 4) imargin(0 0 1 0) graphregion(margin(l=10 r=10)) graphregion(fcolor(white))  saving(pic2, replace)  

gr combine pic1.gph pic2.gph, ycommon col(2) plotregion(fcolor(white)) graphregion(fcolor(white))   title(`"Dependent Variable: Green Tax Levels"', color(black)) saving(Figure3, replace)

